/****************************************************************
Compile: gcc -D_GNU_SOURCE -O3 -o extended extended.c -lm
Usage: ./extended input-1 timestep $RANDOM N

This program is used to calculate the trajectory for a given network
under the extended rules.

Format of the inputfile:
first row            : N (number of the nodes)
second row           : s(1) s(2) ... s(n) (initial state)
the following N rows : network topology (for example, (i,j)=g means that 
                       node j can activate node i
****************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define THRESHOLD 50000

void print_help ()
{
	printf ("Usage         : ./extended inputfile timestep randomseed N\n");
	printf ("                (./extended input-1 0.1 $RANDOM 100)\n\n");
	printf ("This program will run the dynamics with N different parameters under the extended model and reported the results (convergence or not, if yes, how many flips for each node)\n");
}

/* return 1: convergence */
int dynamics (int N, int network[N][N], double para[N], double timestep, double c0[N], double c1[N], int flips[N])
{
	int i, j, p, state[N], flag;

	for (i = 0; i < N; i++)
	{
		if (c0[i] >= 0.5)
			state[i] = 1;
		else
			state[i] = 0;
	}
	
	flag = 0;

	for (i = 0; i < N; i++)
	{
		p = 0;
		for (j = 0; j < N; j++)
		{
			if (i != j)
				p += state[j] * network[i][j];
			else
				p += network[i][j];
		}

		if (p > 0)
		{
			c1[i] = c0[i] + para[i] * timestep * (1.0 - c0[i]);
			if (c1[i] > 0.99)
				flag ++;
		}
		else if (p < 0)
		{
			c1[i] = c0[i] - para[i] * timestep * c0[i];
			if (c1[i] < 0.01)
				flag ++;
		}
		else
		{
			c1[i] = c0[i];
			flag ++;
		}
	}

	for (i = 0; i < N; i++)
	{
		//printf ("%.2lf ", c0[i]);
		if (c0[i] != c1[j] && (c0[i] - 0.5) * (c1[i] - 0.5) < 0)
			flips[i] ++;
	}
	//printf ("\n");
	
	if (flag == N)
		return 1;
	else
		return 0;

}

int main (int argc, char *argv[])
{
	int i, j, n, N;
	FILE *fp;
	char buf[256];

	if (argc == 1)
	{
		print_help ();
		exit (0);
	}

	fp = fopen (argv[1], "r");
	fgets (buf, 256, fp);
	sscanf (buf, "%d", &N);

	int state[N], network[N][N];

	for (i = 0; i < N; i++)
	{
		fscanf (fp, "%d", &state[i]);
	}
	fgets (buf, 256, fp);

	for (i = 0; i < N; i++)
	{
		fgets (buf, 256, fp);
		for (j = 0; j < N; j++)
		{
			if (buf[2 * j] == 'n')
				network[i][j] = 0;
			else if (buf[2 * j] == 'g')
				network[i][j] = 2;
			else if (i == j)
				network[i][j] = -1;
			else
				network[i][j] = -N * 2;
		}
	}

	fclose (fp);

	double timestep, c0[N], c1[N], para[N];
	int flips[N];
	long int seed;

	timestep = atof (argv[2]);
	seed = atoi (argv[3]);
	srand48 (seed);

	for (n = 0; n < atoi (argv[4]); n++)
	{
		for (i = 0; i < N; i++)
		{
			c0[i] = state[i];
			flips[i] = 0;
			para[i] = exp10 (drand48 () - 0.5);
			printf ("%.2lf ", para[i]);
		}
        
		for (i = 0; i < THRESHOLD; i++)
		{
			if (dynamics (N, network, para, timestep, c0, c1, flips) == 1)
				break;
        
			for (j = 0; j < N; j++)
			{
				c0[j] = c1[j];
			}
		}
        
		if (i == THRESHOLD)
			printf ("is not convergenced\n");
		else
		{
			for (i = 0; i < N; i++)
				printf ("%d ", flips[i]);
			printf ("\n");
		}
	}
}
